Quantum Simulation and Phase Diagram of the Transverse Field Ising Model with Three Atomic 

Spins 
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We perform a quantum simulation of the Ising model with a transverse field using a collection of three 
trapped atomic ion spins. By adiabatically manipulating the Hamiltonian, we directly probe the ground state 
for a wide range of fields and form of the Ising couplings, leading to a phase diagram of magnetic order in 
this microscopic system. The technique is scalable to much larger numbers of trapped ion spins, where phase 
transitions approaching the thermodynamic limit can be studied in cases where theory becomes intractable. 
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two energy spectra as a function of the scaled transverse field 
By/\Ji\ in the case of ferromagnetic (FM) nearest-neighbor 
interactions ( Ji < 0). We prepare the system in the ground 
state of the transverse field (By ^^ \Ji\), depicted by the solid 
circle (Fig. [T]), and then adiabatically lower the field com- 
pared to the Ising couplings. When By/\Ji\ <Cl the Ising 
interactions determine both the form of the ground state and 
the energy spacing Age to the excited state(s). 

In Fig. [T^, the next-nearest-neighbor interaction is also FM 
( J2 < 0). There are no level crossings with the ground state 
over the trajectory indicated by the arrow, thus if the evolution 
of the Hamiltonian is slow enough, the system remains in the 
ground state (solid black line). We can also change the sign 
of By and adiabatically follow the highest excited level, as it 
exhibits the same structure as the ground state. In Fig. [TJd, 
J2 > and the next-nearest-neighbor interaction is antiferro- 
magnetic (AFM). The gap at the crossover to magnetic order 
defined by the Ising couplings is ~15 times smaller than that 
of Fig. [It, requiring a slower change of By/ \Ji\ to remain in 
the ground state. 

The competition between different parameters in Eq. [T] 
gives rise to a complex phase diagram. The 2^ possible spin 
configurations are defined as two FM states, |ttt) ^nd |m), 
two symmetric AFM states, |iti) ^^^ I tit)' ^^d four asym- 
metric AFM states, |ttt). Ittt). Ittt) and |||t). all defined 
along the x-axis of the Bloch sphere. In Fig. [2j we plot a part 
of the theoretical phase diagram where the nearest-neighbor 
interactions are always FM ( Ji < 0). The order parameter is 
the probability of occupying an FM state, P(FM) = ^|ttt> + 
P\lll)- For regions where By/\Ji\ ^1, the ground state 
is polarized along By with P(FM)= 1/2^-^ = 1/4. As 
By/\Ji\ decreases, different magnetic phases arise. When the 
next-nearest-neighbor interaction is also FM ( J2 < 0), and 
By/\Ji\ <C 1 the ground states are the two degenerate FM 
states (Fig. [T^). In the region where the next-nearest-neighbor 
interaction is AFM (J2/I Ji| > 0) and J2 overpowers Ji, the 
asymmetric AFM states are lowest in energy. A special point 
next-nearest-neighbor interaction, with aa^ the Pauli spin op- appears at J2/I Ji| = —1 and By = 0, where all the contours 
erator of the i-th particle along the a-direction. Fig. [T] shows of constant FM order meet. Here, the ground state will be a 



At the pinnacle of quantum information science is the 
full scale quantum computer KU, where applications such 
as Shor's factoring algorithm [2J can provide an exponential 
speedup compared with any known classical approach. While 
large-scale quantum computers may not be available for some 
time 1 3 1, more restricted quantum computers known as quan- 
tum simulators look promising right now |4|. As first consid- 
ered by Richard Feynman |[5i|, a quantum simulator controls 
interacting quantum bits (qubits) to implement evolution ac- 
cording to a known Hamiltonian I6j|. Then, by performing 
particular correlation measurements on the qubits, properties 
of certain Hamiltonians - like their ground state - can be ex- 
tracted, often more efficiently than any classical simulation of 
the underlying quantum system |7|. A good example is a col- 
lection of interacting magnetic spins, where the Hamiltonian 
can easily be written down, yet the ground state of magnetic 
order cannot always be predicted, even with just a few dozen 
spins 1 81. 

In this Letter we simulate the Ising model with an trans- 
verse magnetic field and generate a phase diagram, using a 
system of A^ =3 trapped atomic ions. By adiabatically ma- 
nipulating the Hamiltonian, we extract the phases of magnetic 
order in the ground state as a function of the transverse field 
and Ising couplings (SHDj . While this system admits an exact 
theoretical treatment, it also represents the smallest possible 
system having multiple Ising couplings, which give rise to in- 
teresting magnetic order in the phase diagram. Furthermore, 
the experiment is scalable to larger numbers of spins where 
theoretical predictions become intractable. 

The system is described by the Hamiltonian 



H = ^X. 
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with Ising couplings Jij between spins i and j and a uni- 
form transverse magnetic field By. For three spins along a 
symmetric ID chain, we define Ji = Ji 2 = ^2,3 as the 
nearest-neighbor interaction strength and J2 = Ji 3 as the 



superposition of the FM and asymmetric AFM states. This 
effect arises because the pairwise interaction energy cannot 
be minimized individually, leading to a highly degenerate, or 
frustrated, ground state LIU . 
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Figure 1. Energy level diagrams for Eq. [TJwith two different types of 
spin-spin interactions. For both panels, the nearest-neighbor interac- 
tions are FM ( Ji < 0). (a) The next-nearest-neighbor interaction is 
FM with J2/I Ji| ^-2 and (b) AFM with J2/I Ji| --0.9. The arrow 
in both diagrams indicates the trajectory in the simulation, initialized 
at By/\Ji\ ^10. Under this condition, the initial ground state is an 
eigenstate of second term in Eq. [T], a polarized state along By. In 
both examples, at B ^ | Ji | some high energy states cross, but the 
ground state (black solid line) has no level crossings with any excited 
state. Likewise, the highest energy state does not cross any other lev- 
els, allowing one to also adiabatically follow this state. The dotted 
lines represent excited states which are coupled to the ground state 
along the path. In the large field limit, the energy difference between 
ground and excited states Age (here, scaled by y^B^ + Ji) is pro- 
portional to By, but as By/\Ji\ decreases the spin-spin couplings 
determine the energy difference and the form of the ground state. In 
both (a) and (b), the final ground state is FM (defined along the x- 
axis of the Bloch sphere), however in the case of (a), the gap to the 
nearest allowed excited state at the crossover is ^^15 times larger. 

We experimentally simulate the transverse field Ising model 
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Figure 2. Theoretical phase diagram for Eq. [T] The color scale 
indicates the amplitude of the FM order parameter, P{FM) = 
^Ittt) + ^IW- Here, Ji is always negative, yielding FM order in 
that coupling. In the region where J2/I Ji | <0, there is a crossover 
to FM order as By/\Ji\ is lowered (corresponds to Fig. fit). When 
J2/\Ji\ >0, the AFM and FM interactions compete (also shown in 
Fig. [It). This gives rise to a special sharpened point at J2/I Ji | =-1 
and By =0. Here, the ground state is comprised of 6 states: four 
asymmetric AFM and two FM states. 



near the ground state and deeply within the Lamb-Dicke limit 
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The effective magnetic field By, which produces Rabi os- 
cillations between the two spin states, is generated by uni- 
formly illuminating the ion chain with two Raman laser 
beams having a difference frequency at the hyperfine splitting, 
i^HF = 12.642821 GHz. For an individual beam detuning of 
~1.8 THz below the'^Si/2 —'^ P1/2 transitioning and a peak 
intensity of 10 W/ mm^ each ion undergoes Rabi oscillations 
at a rate of 1] ~ 1 MHz and experiences a ~ 20 kHz differ- 
ential AC Stark shift. 

The spin- spin interaction Jij is created by coupling the 
ions' spin states through the normal modes of motion of the 
chain. The two Raman beams travel perpendicular to each 
other to have a wavevector difference Sk along the trans- 
verse direction. The laser frequency of one of the two path- 
ways is modulated to yield beatnotes (with respect to the non- 
modulated beam) at frequencies vhf i M, imparting a spin- 
dependent force at frequency /i (T5\\l6\ . By controlling the 
beatnote detuning /i, we tailor the Ising couplings according 
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of Eq. 1 using cold trapped ^^^Yb^ ions lfT2ll . with the ef- Here, Qi is the Rabi frequency of the i^^ ion. The 



fective spin- 1/2 system represented by the hyperfine "clock" 
states '^Si/2 \F = 1, ttif = 0) and \F = 0, ttif = 0) denoted 
as It)^ and ||)^, respectively (131, where |t)^=|t) + ||) and 
14^)^= It) — It). We confine A^ = 3 atomic spins in a lin- 
ear radiofrequency ion trap and couple them through N col- 
lective transverse motional normal modes along one princi- 
pal axis. These vibrational normal modes, having frequencies 
(^1, i^2, ^3) ^ (4.334, 4.074, 3.674) MHz, are each cooled to 



Lamb-Dicke parameter for the m^^ mode of the i^^ ion is 

Vi,m = bi^rn^ky/h/{4:7rMh'rn), whcrc b IS the normal mode 
transformation matrix and M the mass of a single ion. In 

— /i| >> r]Qi, so that 
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the above expression, we assume 
phonons are only virtually excited. 

We initialize the spins along the 5^ -direction through opti- 
cal pumping (~ 1 /is) and a 7r/2 rotation about the — x-axis of 
the Bloch sphere. The simulation begins with a simultaneous 



and sudden application of both By and Jij where By over- 
powers Jij (By/\Jij\ ^ 10). A typical experimental ramp 
of By decays as By = ae~^^^ + h with a time constant of 
r ~ 30 /iS, varying from a ~ 10 kHz to a final offset of 
h r^ 500 Hz after t = 300 /iS. By varying the power in only 
one of the Raman beams, this procedure introduces a change 
in the differential AC Stark shift of less than 2 Hz. We turn off 
the Ising interactions and transverse field at different By / \ Jij 
endpoints along the ramp. We then measure the magnetic or- 
der along the x-axis of the Bloch sphere by first rotating the 
spins by 7r/2 about the y-3xis, and detecting the 2;-component 
of the spins through spin-dependent fluorescence |[T3l . By re- 
peating identical experiments ~ 1000 times, we obtain the 
probability for the system to be in a particular spin configu- 
ration. We collect fluorescence with a photomultiplier tube, 
exhibiting ~ 97% detection effeciency per spin after 0.8 ms 
exposure. 

This procedure is performed for nine different combina- 
tions of Ji and J2 set by the beatnote detuning /i from Eq. [2] 
In Fig. [3] we present the results as a 3D plot of the FM or- 
der parameter, with the theoretical phase diagram (surface) in 
Fig. [2] superimposed on the data. The data is in good agree- 
ment with the theory (average deviation per trace is ~ 0.09) 
and shows many of the essential features of the phase diagram. 
As By/ \Ji\ decreases, a smooth crossover from a non-ordered 
state to FM order occurs in the region where J2 / 1 Ji | < (Fig. 
[Sj)). As the number of spins increases, this is an example of 
a quantum phase transition. A first order transition due to an 
energy level crossing is apparent (Fig. [3]:) when changing J2 
for a fixed and small value of By/ \Ji\ = 0.57. This tran- 
sition is sharp, even in the case of three spins. The data (e.g. 
Fig. [Sj)) show small amplitude oscillations in the initial evolu- 
tion due to the sudden application of the spin- spin interaction, 
which is held constant during the simulation to minimize vari- 
ation in the differential AC stark shifts. This limitation can be 
removed by choosing the Raman laser detuning such that con- 
tributions from the ^Ps/2 energy level lead to a minimum in 
the ratio of the differential AC stark shift to the resonant Rabi 
frequency ifTTl . 

We now investigate adiabaticity of the Hamiltonian trajec- 



tory H{t), characterized by the condition iflS I 
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In this expression. 



the dimensionless quantity 

characterizes coupling from the 

ground state \g{t)) to any excited state \e{t)) with energy 
gap hAge. This parameter is small, of order unity for this 
simulation, but is peaked at a crossover in magnetic order, 
where the instantaneous eigenstates are most rapidly varying. 
Therefore, Eq. |3] states that to remain adiabatic, the slope of 
the time-dependent By -fiQld profile must be shallow when 
the gaps in the energy spectrum are small (as in Fig. [TJ3), in 
particular near a crossover (phase transition for large N). 
In Fig. |4j we investigate this adiabatic criteria for two dif- 
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Figure 3. Experimental measurements of the phase diagram for 
Eq. [T](solid bars) compared to the theoretical prediction from Fig. |2] 
(surface) . The vertical amplitude is the FM order parameter P(FM)= 
P\-^-^-^) + P III) . The ratio of By/\Ji \ was varied from ^10 to ^0.1 
for J2/I Ji values of -1.3,-2.0, -3.6, 4.2, 2.0, 1.3, 0.92, 0.74, and 
0.62. Ji < for all traces, (b) As By/ \Ji\ ^ in the region 
where J2/IJ1I < — 1 , we observe a smooth crossover to FM or- 
der. The filled circles and solid line are the data and theory for 
J2/\Ji\ = —1.3, respectively (c) When changing J2 for a fixed 
and small value of By/ \Ji\ the system undergoes a sharp transition. 
The data (filled circles) shown is for a scan of By/ Jy = 0.57. The 
average deviation per scan of By/ \Ji\ from the exact ground state 
is - 0.09. 



ferent types of next-nearest-neighbor coupling. In Fig. |4^ all 

^ -2 (as in Fig. [T^). The 



interactions are FM and J2/I Ji 
dashed lines in the top panel are the adiabaticity parameter 
from Eq. |3] calculated over the trajectory for the two coupled 
excited states (recall Fig. [TJ. Due to the 500 Hz final offset of 



By, the simulation stops atBy/ Ji '^ 0.5. To examine the be- 
havior extended below this value, we calculate the criteria for 
an exponential ramp with a 100 /is time constant. This pro- 
file was chosen to overlap with experimental parameters for 
large By/\Ji\ and also reach By = in a typical simulation 
time (~ 300 /is). The results indicate that Eq. 3] is satisfied 
over the trajectory; By{t)e/A'^^ remains much less than one 
even with a maximum occurring Sit By/ \Ji\ ~ 1. To demon- 
strate the simulation is indeed adiabatic for these parameters, 
the measured probability P(FM) (solid dots) is shown in the 
lower panel of Fig. |4^. The black line represents the adia- 



batic ground state and the grey line is the theoretical expected 
probability including the experimental ramp. The dotted line 
in this figure is the theoretical state evolution using a By -fiQld 
ramp that reaches zero. The predicted evolution does not sig- 
nificantly deviate from the ideal ground state and the data is 
in good agreement with all three theory curves. 

Fig. |4]3 presents the case when the next-nearest-neighbor 
interaction is AFM and J2/I Ji| ~ 0.9 (as in FiglTJ)). When 
By/\Ji\ <C 1, By{t)e/ /S?g^ reaches a maximum value of 
~ 0.6, indicating that the probability for excitations will likely 
increase. This difference is because in this case the gap A^e at 
the 'critical' point is ~ 15 times smaller than that in Fig. [T^. In 
contrast to the FM J2 case, the theoretical probability curves 
shown in the lower panel of Fig. [43 predict significant dia- 
batic effects when using this 5^ -field profile for simulations 
near the special point. In fact, to successfully evolve to the 
true ground state near By = 0, the simulation time (assuming 
same initial conditions and an exponential ramp of By) should 
be at least a factor of ten longer. 

Because all the data lies outside of the region where the en- 
ergy gaps are small, the diabatic excitations are minimal, but 
further experimental study is needed to precisely quantify this 
effect. One method to probe excitations, which may also be 
useful as A^ ^ 1, is to perform and then reverse the experi- 
mental ramp and measure the probability of returning to the 
initial state. The main contributions to the overall data off- 
set from theory shown in Fig. [4] and Fig. |3] are spontaneous 
emission due to off-resonant scattering (probability ~ 15% in 
1 ms), imperfect optical pumping (state preparation), parasitic 
fields along the x— and z— axes, and state detection error. Ad- 
ditional phonon terms not appearing in the Hamiltonian of Eq. 
[T]are expected to contribute at a level under 2% ifTOl . 

As the number of spins N grow, the technical demands on 
the apparatus are not forbidding |10|[TT1. In particular, the 
expected adiabatic simulation time for this model is inversely 
proportional to the 'critical' gap in the energy spectrum; for 
the transverse field Ising model in a finite- size system, this gap 
decreases as A""^/^ (19J. Scaling this system to accommo- 
date long ion chains (approaching the thermodynamic limit) 
will allow investigation of behavior near critical points. This 
is interesting for A' >20, where general spin models become 
theoretically intractable. For instance, the Lanczos algorithm 
(20 1 can be used to find low-lying states of a 30 spin/site prob- 
lem if the 2^^ element matrix is sparse. If one is interested in 
a non- sparse matrix, as is the case for the ion system with 
long-range magnetic coupling, the limiting number of spins is 
^20. Theoretical investigations of dynamics limit the number 



of spins further ||81[2TJ. We note that this approach to quan- 
tum simulation is versatile and may be extended to simulate 
Heisenberg or XYZ spin models using additional laser beams 
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Figure 4. Adiabaticity for two cuts in the phase diagram where (a) 
J2/\Ji\ = —2; and (b) J2/IJ1I = 0.9. The upper plots show the 
theoretical ground state adiabaticity parameter By (t)e/ A^^ of Eq . |3| 
(dashed lines) for each of the two coupled excited states (see Fig. [ij. 
We use a pure exponential decay ramp for the field By with time 
constant 100 /is in order to match the experiment for large values of 
By/Ji while also extending the theory curves below the minimum 
value of By reached in the experiment. At every point in the the 
trajectory of (a), we find that the adiabatic condition (Eq. |3]) is satis- 
fied for typical experimental times (300 /xs). On the other hand, in 
(b), there is a significant probability of diabatic transitions to excited 
states for By/Ji <^ 1. In the lower plots, we compare the observed 
FM order parameter (points) with theory. In (a), the theoretical order 
from the exact experimental ramp with a 35 /xs time constant and 
final offset value given in the text (grey solid line) is in reasonable 
agreement with the order in the true ground state (black solid line) 
for By/ Ji > 0.5. The dotted line is the expected state evolution for 
a pure exponential decay ramp with a 100 /us time constant, allowing 
By ^ 0. In (b), the data also matches well to theory, as we avoided 
the regions where diabatic transitions are expected for By/Ji ^ 1. 
According the calculations, the duration of three-spin experiments 
near the special point should be on the order of milliseconds. 
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